/*-----------------------------------------------------------------------------
 * File Name ..... util.c
 *
 * Purpose ....... Contains functions used for gridding
 *
 * Functions ..... MyMalloc, MyCalloc, mysecond, Multiply
 *
 * Date Written .. Nov 30, 2010
 *
 *---------------------------------------------------------------------------*/

#define _GNU_SOURCE
#include <sched.h>
#include <stdio.h>
#include <sys/syscall.h>
#include <sys/time.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <pthread.h>
#include "config.h"
#include "userdef.h"

#ifndef PAGESIZE
#define PAGESIZE    ( 4 << 10 )
#endif
#define HUGE_PAGESIZE   ( 2 << 20 )

#define CACHE_PROP_CS 1.0
#define CACHE_PROP_CD 0.7


#ifdef ROWMAJOR
#define A(I,J) A[((I))*n+(J)]
#define B(I,J) B[((I))*n+(J)]
#define C(I,J) C[((I))*n+(J)]
#define D(I,J) D[((I))*n+(J)]
#else
#define A(I,J) A[((J))*n+(I)]
#define B(I,J) B[((J))*n+(I)]
#define C(I,J) C[((J))*n+(I)]
#define D(I,J) d[((J))*n+(I)]
#endif



#define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
#ifndef max
#define max( a, b ) ( ((a) > (b)) ? (a) : (b) )
#endif

/*----------------------Functions--------------------------------------------*/
void *MyMalloc(size_t MemSize);
void *MyCalloc(size_t MemSize);
void PrintData(double *Array, size_t DimensionCount);
double get_time();
void *Multiply(void *WorkerData);
extern pthread_mutex_t GridMutex;	/* Mutex for grid updation*/
void fill(double *Matrix, int size);
/*--------------------------------------------------------------------------*/

/*----------------------external Variables--------------------------------------------*/
extern double TotalTime;
extern int ThreadCount; /* No of Threads */
extern pthread_barrier_t Barrier;
/*--------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------
* Function ........... MyMalloc
*
* Purpose ............ Allocates memory using MyMalloc and does error checking
*
* Arguments Passed ... MemSize - size of memory to be allocated
*
* Return Value ....... Pointer to the allocated memory
*
* Globals Modified ... None
*---------------------------------------------------------------------------*/
void *MyMalloc(size_t MemSize)
{
    void *ptr;

    if ( (ptr = (void *) malloc(MemSize)) == (void *) NULL)
    {
        fprintf(stderr, "Error (MyMalloc) : Unable to allocate memory\n");
        exit(0);
    }
    return ptr;
}


/*----------------------------------------------------------------------------
 * Function ........... MyCalloc
 *
 * Purpose ............ Allocates memory using calloc (for 1 element only)
 *                      and does error checking
 *
 * Arguments Passed ... MemSize - size of memory to be allocated
 *
 * Return Value ....... Pointer to the allocated memory
 *
 * Globals Modified ... None
 *---------------------------------------------------------------------------*/
void *MyCalloc(size_t MemSize)
{
    void *ptr;

    if ( (ptr = (void *) calloc(1, MemSize)) == (void *) NULL)
    {
        fprintf(stderr, "Error (MyCalloc) : Unable to allocate memory\n");
        exit(0);
    }
    return ptr;
}



/*----------------------------------------------------------------------------
 * Function ........... PrintData
 *
 * Purpose ............ prints the values stored in the arrays
 *
 * Arguments Passed ... array size and a pointer to array's base
 *
 * Return Value ....... none
 *
 * Globals Modified ... None
 *---------------------------------------------------------------------------*/
void PrintData(double *Array, size_t size )
{
    size_t item = 0;

    while (item < size)
    {
        fprintf(stdout, " %10.10f\n", Array[item++]);
    }
    fprintf(stdout,"\n------------------------------------\n");
}


/*----------------------------------------------------------------------------
* Function ........... get_time
*
* Purpose ............ Returns the time (seconds+microseconds) based on gettimeofday
*
* Arguments Passed ... none
*
* Return Value ....... double that has the time
*
* Globals Modified ... None
*---------------------------------------------------------------------------*/
double get_time()
{
    struct timeval tp;
    struct timezone tzp;
    int i;
    i = gettimeofday(&tp,&tzp);
    return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
}

void fill(double *Matrix, int size)
{
    int i;
    for (i=0; i<size; i++)
    {
        Matrix[i] = rand() % 2000;
    }
}

void my_dgemm(const int* m,const int* n,const int* p, const int *alpha, const double *a,const int* lda,const double * b,const int* ldb,const int* beta,double * c,const int* ldc) {
    int i,j,k;
    for (i=0; i<*m; ++i) {
        for (j=0; j<*n; ++j) {
            double tmp = c[i**ldc+j];
            for (k=0; k<*p; ++k) {
                tmp += a[i**lda + k]* b[k**ldb + j];
            }
            c[i**ldc+j] = tmp;
        }
    }
};

#define dgemm(CblasATrans,CblasBTrans,m,n,p,alpha,a,lda,b,ldb,beta,c,ldc) my_dgemm(m,n,p,alpha,a,lda,b,ldb,beta,c,ldc)

#define cblas_dgemm(CblasXMajor,CblasATrans,CblasBTrans,m,n,p,alpha,a,lda,b,ldb,beta,c,ldc) my_dgemm(&m,&n,&p,&alpha,a,&lda,b,&ldb,&beta,c,&ldc)




int get_lambda(int cs) {
    return max((int) floor(sqrt((double)cs -0.75) - 0.5),1);
}

int get_alpha_max(int cs) {
    return max((int) floor(-1.0 + sqrt((double)cs + 1.0)),1);
}


int get_mu(int cd) {
    return max((int) floor(sqrt((double)cd -0.75) - 0.5),1);

}

#define ADD_INV_MOD_P(x,p) (((int)(p) -( (int)(x) % (int)(p))) % (int)(p) )
void get_alpha_beta(int cs, double bws, double bwd, int p,int mu, int block,int  n,int * alpha, int * beta) {
    *alpha = floor(sqrt((double)cs*(1.0+2.0*(double)p*bwd/bws - sqrt(1.0+8.0*(double)p*bwd/bws))/(2.0*(double)p*bwd/bws -2.0)));

    *alpha += ADD_INV_MOD_P(*alpha,(int)sqrt(p)*mu);
    *alpha -= (int)sqrt(p)*mu;

    while (*alpha>sqrt(p)*mu && (*alpha%((int)sqrt(p)*mu)!=0 || n%(*alpha*block)!=0) )
    {
        --*alpha;
    }


    *beta = ((double)cs - (double)pow(*alpha,2))/(4.0**alpha);
    while (*beta>1 && n%(*beta*block)!=0  )
        --*beta;
}




/*----------------------------------------------------------------------------
* Function ........... Multiply
*
* Purpose ............ Imposes a grid on k-space data. Gridding is done according to the formula
*                      gridpoint value += datavalue/dcf * kernel(dk)
*
*
* Arguments Passed ... x and y values of samples, real and imaginary values at x,y
*		        , Number of samples, dcf array, pointer to output arrays (real and
* 			imag), gridsize, kernel width in grid points, conv kernel values from
*			0 to its width, number of points in kernel lookup table
*
* Return Value ....... none
*
*
* Globals Modified ... output arrays
*---------------------------------------------------------------------------*/
void *Multiply(void *WorkerData)
{
    THREAD_DATA *my_data; 			/* To store each threads data */
    const DATA_TYPE *A, *B;
    DATA_TYPE *C;
    int A_Rows, A_Cols, B_Rows, B_Cols, C_Rows, C_Cols;
    int ThreadId; 				/* Thread's Id to find it's own data range */
    int Size_A = 0, Size_B = 0, Size_C = 0;
    int steps;
    int n;
    int mu;
    int sqrtproc;
    int tmp;
    int blk_size;
    int alphamax,alphamin,alpha,beta;
    int alphap, maxstep,maxsubij;
    int step,maxij;
    int i,j,p,stepi,stepj,k,ia,substep;
    int my_i,my_j,subi,subj,ka,ja,hdim;
    int numproc, mn, blocksize;
    int cs, cd;

    int it,jt,ig,jg,im,jm,ib,jb;

    /* Read arguments into local pointer */
    my_data = (THREAD_DATA *) WorkerData;

    /* Thread should know its ID */
    ThreadId = my_data->ThreadId;

    /* Read all args from arguments */
    A      = my_data->A;
    B      = my_data->B;
    C      = my_data->C;
    A_Rows = my_data->A_Rows;
    B_Rows = my_data->B_Rows;
    C_Rows = A_Rows;
    A_Cols = my_data->A_Cols;
    cs     = my_data->cs;
    cd     = my_data->cd;
    A_Cols = my_data->A_Cols;
    B_Cols = my_data->B_Cols;
    C_Cols = B_Cols;
    blocksize = my_data->blocksize;
    numproc = my_data->numproc;
    double coms = 2429340876;
    double comd = 4563402752 ;


    /*copying parameters in private variables*/

    mn       = A_Rows * B_Rows;
    n        = mn;


    //need to find parameters lambda and mu
    mu = get_mu(cd);

    sqrtproc = (int)floor(sqrt(numproc));
    tmp = mu;
    while (tmp>0 && n%(tmp*sqrtproc*blocksize) !=0 )
    {
        tmp--;
    }

    mu = max(tmp,1);

    alphamax = get_alpha_max(cs);
    alphamin = sqrtproc*mu;
    get_alpha_beta(cs, coms, comd, numproc, mu, blocksize,n,&alpha,&beta);

    alphap = alpha/sqrtproc;

    int alphablk = alpha*blk_size;
    int giga_size = alpha/(sqrtproc*mu);
    int mega_size = mu;


    maxsubij = max(alpha/alphamin,1.0);

    step = alpha;
    hdim = mu*blk_size;

    maxij = (int)ceil((double)n/((double)step*(double)blocksize));
    maxstep = (int)ceil((double)n/((double)beta*(double)blocksize));
    #ifdef ROWMAJOR
    my_i =  (ThreadId % sqrtproc);
    my_j =  (ThreadId / sqrtproc);

    for (stepi = 0; stepi <maxij; ++stepi)
    {
        i = stepi*step;
        for (stepj = 0; stepj <maxij; ++stepj)
        {
            j = stepj*step;
            for (substep=0; substep<maxstep; ++substep) {

                k= substep*beta;

                for (subi = 0; subi <maxsubij; ++subi)
                {
                    for (subj = 0; subj <maxsubij; ++subj)
                    {
                        for (ka = k; ka <k + beta; ++ka)
                        {
                            for (ia = i + my_i*alphap + subi*mu ; ia<i + my_i*alphap + (subi+1)*mu; ++ia)
                            {
                                double one = 1.0;
                                for (ja = j + my_j*alphap + subj*mu ; ja<j + my_j*alphap + (subj+1)*mu; ++ja)
                                {
                                    cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,blk_size,blk_size,blk_size,one,&A(ia*blk_size,ka*blk_size),n,&B(ka*blk_size,ja*blk_size),n,one,&C(ia*blk_size,ja*blk_size),n);
                                }
                            }
                        }
                    }
                }
            }
        }
    }



#else
    my_i =  (ThreadId % sqrtproc);
    my_j =  (ThreadId / sqrtproc);

    for (stepj = 0; stepj <maxij; ++stepj)
    {
        j = stepj*step;

        for (stepi = 0; stepi <maxij; ++stepi)
        {
            i = stepi*step;

            for (substep=0; substep<maxstep; ++substep) {

                k= substep*beta;
                for (subj = 0; subj <maxsubij; ++subj)
                {

                    for (subi = 0; subi <maxsubij; ++subi)
                    {
                        for (ka = k; ka <k + beta; ++ka)
                        {
                            for (ia = i + my_i*alphap + subi*mu ; ia<i + my_i*alphap + (subi+1)*mu; ++ia)
                            {
                                for (ja = j + my_j*alphap + subj*mu ; ja<j + my_j*alphap + (subj+1)*mu; ++ja)
                                {
                                    dgemm("N","N",&blk_size,&blk_size,&blk_size,&alpha,&A(ia*blk_size,ka*blk_size),&n,&B(ka*blk_size,ja*blk_size),&n,&alpha,&C(ia*blk_size,ja*blk_size),&n);
                                }

                            }
                        }
                    }
                }
            }
        }
    }
#endif

    pthread_exit(NULL);

}

